!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate inverse wall 
! papramenters by MOC.
! ref. to chapter 15 of 'Gas Dynamics' by Maurice J. Zucrow & Joe D. Hoffmann
! -----Sauer Method------
! -----Revised by Kliegel Method
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-06-24
! revised  : 2015-06-26    Kliegel Method for accuracy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : Transonic
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   epsi      displacement of coordinate origin relative to the nozzle throat
!   astar     sonic speed of stagnation condition
!   alpc      coefficient of linear axial perturbation velocity
!   C         coefficients for calculate mdot and F
!   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   1. angle in rad, not degree

subroutine Transonic(xbarma,xbarv,y,Vuwxy,Vvwxy,V,M,T,p,rou)
   use VariableDef
   implicit none
!     given variable
!   real*8::
!     temp variable
   real*8::h,epsi,astar,alpc,roustar,pstar,Tstar,xma(Ny),xv(Ny),C(Ny)
   real*8::Rcap,zz,rr,ut1,ut2,ut3,vt1,vt2,vt3
   integer::i
!     calculate variable
!   real*8::mdot,F,mdot1d,F1d,CDm,lambdaF
   real*8::y(Ny),xbarma(Ny),xbarv(Ny),Vuwxy(Ny),Vvwxy(Ny),V(Ny),M(Ny),T(Ny),p(Ny),rou(Ny)

   epsi=-yt*sqrt((g+1)*(1+delta)/(rtu/yt))/(2*(3+delta))
   astar=sqrt(2*g*Rg*Tt/(g+1))
   Tstar=2*Tt/(g+1)
   pstar=Pt*(2/(g+1))**(g/(g-1))
   roustar=pstar/(Rg*Tstar)
   alpc=sqrt((1+delta)/((g+1)*rtu*yt))
   Rcap=rtu/yt

   h=yt/(Ny-1)

!    calculate properties along the initial line   
   do i=1,Ny
      y(i)=h*(i-1)   ! throat height step
      xma(i)=-alpc*(g+1)*y(i)**2/(2*(1+delta))   ! Mach 1 line axial location
      xv(i)=-alpc*(g+1)*y(i)**2/(2*(3+delta))   ! zero radical velocity line axial location
      xbarma(i)=xma(i)-epsi   ! Mach 1 line axial location at throat coordinate system
      xbarv(i)=xv(i)-epsi   ! zero radical velocity line axial location at throat coordinate system

!      corrected by Kliegel method
      rr=y(i)/yt
      zz=sqrt(2*Rcap/(g+1))*xbarv(i)/yt
      
      ut1=0.5*rr**2-0.25+zz
      ut2=(2*g+9)*rr**4/24-(4*g+15)*rr**2/24+(10*g+57)/288+zz*(rr**2-5/8)-(2*g-3)*zz**2/6
      ut3=(556*g**2+1737*g+3069)*rr**6/10368-(388*g**2+1161*g+1881)*rr**4/2304+(304*g**2+831*g+1242)&
         &*rr**2/1728-(2708*g**2+7839*g+14211)/82944+zz*((52*g**2+51*g+327)*rr**4/384-(52*g**2+75*g&
         &+279)*rr**2/192+(92*g**2+180*g+639)/1152)+(-(7*g-3)*rr**2/8+(13*g-27)/48)*zz**2+&
         &(4*g**2-57*g+27)*zz**3/144
      vt1=0.25*rr**3-0.25*rr+rr*zz
      vt2=(g+3)*rr**5/9-(20*g+63)*rr**3/96+(28*g+93)*rr/288+zz*((2*g+9)*rr**3/6-(4*g+15)*rr/12)&
         &+rr*zz**2
      vt3=(6836*g**2+23031*g+30627)*rr**7/82944-(3380*g**2+11391*g+15291)*rr**5/13824&
         &+(3424*g**2+11271*g+15228)*rr**3/13824-(7100*g**2+22311*g+30249)*rr/82944&
         &+zz*((556*g**2+1737*g+3069)*rr**5/1728-(388*g**2+1161*g+1881)*rr**3/576&
         &+(304*g**2+831*g+1242)*rr/864)+((52*g**2+51*g+327)*rr**3/192&
         &-(52*g**2+75*g+279)*rr/192)*zz**2-((7*g-3)*rr/12)*zz**3
      
      Vuwxy(i)=astar*(1+ut1/(Rcap+1)+(ut1+ut2)/((Rcap+1)**2)+(ut1+2*ut2+ut3)/((Rcap+1)**3))
      Vvwxy(i)=astar*sqrt((g+1)/(2*(Rcap+1)))*(vt1/(Rcap+1)+(3*vt1/2+vt2)/((Rcap+1)**2)&
         &+(15*vt1/8+5*vt2/2+vt3)/((Rcap+1)**3))
      V(i)=sqrt(Vuwxy(i)**2+Vvwxy(i)**2)
      M(i)=V(i)/sqrt(g*Rg*Tt-(g-1)*V(i)**2/2)

!      end corrected by Kliegel method
      
!      Vuwxy(i)=astar*(1+alpc*xv(i)+(g+1)*alpc**2*y(i)**2/(2*(1+delta)))
!      Vvwxy(i)=0.0
!      V(i)=sqrt(Vuwxy(i)**2+Vvwxy(i)**2)
!      M(i)=V(i)/sqrt(g*Rg*Tt-(g-1)*V(i)**2/2)
      T(i)=Tt/(1+(g-1)*M(i)**2/2)
      p(i)=Pt/((1+(g-1)*M(i)**2/2)**(g/(g-1)))
      rou(i)=p(i)/(Rg*T(i))
   end do
   
!    calculate coefficient Ci
   C(1)=1.0
   do i=2,Ny-1
      C(i)=3.0+(-1)**i
   end do
   C(Ny)=1.0

!    calculate mass flow rate and thrust cross vwbar line
   mdot=0.0
   F=0.0
   do i=1,Ny
      mdot=mdot+2*pi*(h/3)*C(i)*rou(i)*Vuwxy(i)*y(i)
      F=F+2*pi*(h/3)*C(i)*(p(i)+rou(i)*Vuwxy(i)**2)*y(i)
   end do

!   calculate one-dimensional mass flow rate and thrust
   mdot1d=roustar*astar*pi*yt**2
   Fs=pstar*pi*yt**2+mdot1d*astar
   F1d=Fs-p0*pi*yt**2
   CDm=mdot/mdot1d
   etaF=F/F1d
end subroutine Transonic
